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Abstract 


This article comments on a frequency estimator which was proposed by 
[6] and shows empirically that it exhibits a much larger mean squared error 
than a well known frequency estimator by [8]. It is demonstrated that by 
using a heuristical adjustment |2] the performance can be greatly improved. 
Furthermore, references to two modern techniques are given, which both nearly 


attain the Cramér-Rao bound for this estimation problem. 


1 Introduction 


The problem of achieving a precise estimation of the fundamental frequency of a 
noisy signal has been researched for a considerable amount of time. It has many 
applications in different areas of engineering, not least in audio processing and mu- 
sicology. Quite a number of algorithms has been published [5, 4] until today. These 
methods often differ strongly regarding their general idea, area of application, pre- 
ciseness and time complexity. This work deals with an empirical comparison of two 
such algorithms for frequency estimation of simple sinudoids with noise. The basic 
tool for this task is the discrete Fourier transform (DFT) and its efficient implemen- 


tation the fast Fourier transform (FFT). 


*This work has been supported by the Deutsche Forschungsgemeinschaft, Sonderforschungs- 
bereich 475. 


2 Frequency estimation 


Let the signal be a simple sinsusoid 


X,=pt+Asinwot+¢)+e, t=0,1,...,7-1 


with additive noise €,, sampled at T points in time. In general the noise is assumed 
to be Gaussian white noise with Ele] = 0 and Varle] = 0? (although some cited 
papers generalize upon this |7, 8, 9]). While all other parameters are unknown, 
a frequency estimator w is sought, which is as accurate as possible and can be 
computed efficiently in time. Minimizing the squared error seems to be a reasonable 


first approach: 


T-1 
SSK = SOX —p—Asin(wt + 6)? 


t=0 
This equates to a maximum-likelihood approach. The Cramér-Rao bound for unbi- 


ased estimators w is known and amounts to: 


60? 
i ig a SOEs 
Var(w) > TR’ 


which contains apart from the cubic term T(T? — 1) the signal-to-noise ratio (SNR) 
A? /2a? (see [4] for a derivation). 


2.1 Interpolation of Fourier coefficients 


There are quite a number of algorithms for frequency estimation, which more or less 


directly interpolate the complex-valued DFT coefficients 


T-1 


= S- X,exp(—i27jt/T) 


t=0 


or the real-valued coefficients of the periodogram 


T-1 2 


S| X,exp(—itA) 


t=0 
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I(4) = 5 
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where J(A) is in practice only available at the discrete Fourier frequencies 27jt/T: 


_ 2 / 
en | 


The reason for this course of action is that the maximum of the periodogram 


Ie 


= [(2njt/T) 


Wmax = argmax J(\) is an asymptotically unbiased estimator for w, so that also 


asymptotically holds [8]: 


Po) Ges. —w)~N (0 aa , 
where f(w) is the spectral density of the noise e,. However, the global maximum 
of I(A) is hard to compute numerically (in an efficient way) because of the many 
local optima of the oscillating /(\). One can neither use the maximiser of the 
discrete periodogram P, as an approximation - as the resolution of P;, is not good 
enough, even for signals without noise, to get sufficiently precise estimations for 
applications [6] - nor does choosing the maximiser of P, as starting point for an 
iterative optimisation of /(\) guarantee convergence to the global maximum of J()). 
But as the maximum of the periodogram and its two neighbours already contain 
about 85% of the power of the spectrum, an interpolation of these three points 


seems appropriate. 


2.2 Algorithm by Ligges 


Ligges [6] presents a heuristically motivated procedure for frequency estimation in 
the area of musicology. As this algorithm should also be used for signals with 
overtones or multiple tones and implements further practical aspects, it is more 
comprehensive than the following estimators given below, but in the here considered 
case of simple tones it can be reduced to the following: A parabola is fitted to the 
maximum of the periodogram and its neighbouring frequencies and the peak position 


of the parabola is used as an estimator for the tone frequency. 


ee eA OP SOF) IO) a OR. Pt 
age 2. | 2(*) — 10") — 10") TR AP OP A OP io 


Here, A* denotes the maximising frequency of the periodogram, ** the neighbouring 
frequency with the larger and \*** the neighbouring frequency with the smaller value 
in the periodogram. 

Practically the same method was published by Voglewede in [10]. A further heuris- 


tical adaptation to improve the estimation |6, 1] yields: 


LTA) 
2 T(*) 


A * 
WLigges2 = AA 


In his PhD thesis, Ligges [6] conducts an experimental evaluation of these two es- 
timators and reports a mean error of approximately 2.2 Hz for the fundamental 
frequency when using the second algorithm. He also compares his method to known 


model-based approaches. Asymptotic distributions are not indicated. 


2.3 Algorithm by Quinn 


[8] uses an estimator, which in a similar efficient fashion interpolates three Fourier 
coefficients, although he does not use the periodogram but instead directly works 
on the complex DFT coefficients Y; of the data. 


1. Let k be the maximising index of |Y?}. 


2. Let ay = R(Vpp—-1/ Yer) and a2 = R(Vk-41/Ve-), 
and 6; = a;/(1 — a1) and 62 = —a2/(1 — ag) 


3. If both 61,62 > 0, then 6 = do , else 6 = 64 


Quinn does not assume the noise to be i.i.d Gaussian, but derives under quite weak 
and somewhat technical assumptions (e, strictly stationary and ergodic, Ele,| = 0, 
Ele?] < ov, for further details see [8]) the following asymptotic distribution by using 
a central limit theorem: 
73/2 
—— (WQuinn — #) ~ N(O,1) 
UT 


where 


167°6?(1 — |d|)?(26? — 2|5| + 1) f(w) 
A? sin?(16) 
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and f(w) is the spectral density of the noise. 


2.4 Algorithm by Jacobsen 


As the following empirical comparison shows, neither Wrigges1 NOY Wrigges2 Produce 
very good results. Wrigges: can be rectified, if one employs the same general idea as 


in 2.3 and switches to the complex-valued DFT coefficients: 


‘ 4 2m Yp-1 = Yaa 
Jacobsen ?P 2Y, — Yrsi — Ye-i 


Please note the other slight changes in the formula as well. This method was pub- 
lished by Jacobsen in [3]. There are no further theoretical results concerning the 
distribution of the estimator. Also, it’s not very evident, why exactly this adaptation 


results in such an improvement|Jacobsen 2009, personal communication]. 


3 Empirical comparison 


1000 frequencies for w between 2 and 32 generated, while the noise was varied from 
o* = 0 to 1. The amplitude A was set to 1 (so SNR was at least 0.5) and the 
phase @ was selected randomly for the resulting 30000 sinusoids. Every signal was 
sampled 128 times. Figure 1 shows the error distributions for all four estimators. 
It is clearly visible in fig. 1 that the simple quadratic interpolation results in the 
worst accuracy. Both algorithms by Ligges exhibit a much larger variance than the 
methods by Quinn and Jacobsen, also the main mass of the distribution is not even 
centred around zero. Wrigges2 does improve on the variance to some extend, but 
does not change the general, problematic shape of the error function. Jacobsen’s 


estimator seems to have about the same quality as the one by Quinn. 


Fig. 1: Empirical distributions of estimation errors 


To gain some further insights into the estimators, figures 2 and 3 show the devel- 
opment of bias and variance in relation to the true frequency of the signal - first 
without noise, then with o? = 0.2 and at last with 0? = 0.4. While both estimators 
of Ligges are strongly biased and also exhibit a large variance (which of course grows 
with the noise), the remaining two algorithms of Quinn and Jacobsen seem again 
roughly equal in quality, bias and variance being slightly different in different areas 
of the bin. 
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Fig. 2: Bias and variance of the 4 estimators against the bin number of the true 
frequency. The three rows correspond to a noise setting of 0? = 0, 0? = 0.2 and 


ao” = 0.4 respectively. 
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Fig. 3: Bias and variance of best 2 estimators against the bin number of the true 
frequency. The three rows correspond to a noise setting of 0? = 0, 0? = 0.2 and 


o” = 0.4 respectively. 


4 Best known methods near the Cramér-Rao bound 


Apart from the already discussed methods, there a two other ones, which exhibit an 
even lower MSE. The first one is an improvement of Quinn’s algorithm [9] OQuinn, 
where 6; and 6) are combined in a nonlinear fashion to a new 6 and not selected 


depending on their sign: 


_ 01 + 42 
- 2 


1 V6 g+1—,/2/3 
K(x) = —log(3a? + 62 +1) — lo 
(2) d al ) 24 (Hee) 


The second is a further estimator by [7]. MacLeod’s estimator is - on a superfi- 


ft) 


— K(5}) + £(93) 


cial level - similar to Jacobsen’s quadratic estimator plus an additional nonlinear 
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correction term: 
BS hs i; = R(Yi7") ’ 


— Rea Rept 
2Re+ Regi t+ Rear ’ 
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W MacLeod = 
Ay : T 


where k denotes the maximising index of the discrete periodogram and 7* the com- 


oe 


plex conjugate. Though one should note the different signs in the denominator of 
y and Jacobsen’s interpolation quotient and that the Y; are multiplied by a phase 
reference to align them better with the phase of the periodogram peak [7]. 

It can be shown in both cases, that the variance of the estimators is quite close 
to the Cramér-Rao bound. Jacobsen has published his own empirical evaluation of 


these methods on his web page [2]. 


5 Summary 


A couple of methods for frequency estimation of noisy sinusoids - all using the 
interpolation of three Fourier coefficients - were compared. It was demonstrated, 
that both algorithms by Ligges exhibit a large bias and variance. The first estimator 
can be corrected and is then equivalent to an estimator by Jacobsen, which has a 
roughly similar performance as a well known estimator by Quinn. References to two 
further methods were given, which show an even lower MSE than the two mentioned 
before and which are provably nearly optimal with regard to the Cramér-Rao bound. 
In particular, the algorithm by Quinn has the advantage that asymptotic normal 


distributions are available. 
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